N-Heterocyclic Carbene Formation in the Ionic Liquid [EMIM+][OAc–]: Elucidating Solvation Effects with Reactive Molecular Dynamics Simulations

Recent experimental and theoretical work has debated whether N-heterocyclic carbenes (NHCs) are natively present in imidazolium-based ionic liquids (ILs) such as 1-ethyl-3-methylimidazolium acetate ([EMIM+][OAc–]) at room temperature. Because NHCs are powerful catalysts, determining their presence within imidazolium-based ILs is important, but experimental characterization is difficult due to the transient nature of the carbene species. Because the carbene formation reaction involves acid–base neutralization of two ions, ion solvation will largely dominate the reaction free energy and thus must be considered in any quantum chemical investigation of the reaction. To computationally study the NHC formation reaction, we develop physics-based, neural network reactive force fields to enable free energy calculations for the reaction in bulk [EMIM+][OAc–]. Our force field explicitly captures the formation of NHC and acetic acid by deprotonation of a EMIM+ molecule by acetate and in addition describes the dimerization of acetic acid and acetate. Using umbrella sampling, we compute reaction free energy profiles within the bulk IL and at the liquid/vapor interface to understand the influence of the environment on ion solvation and reaction free energies. Compared to reaction of the EMIM+/OAc– dimer in the gas phase, the bulk environment destabilizes formation of the NHC as expected due to the large ion solvation energies. Our simulations reveal a preference for the product acetic acid to share its proton with an acetate in solution and at the interface. We predict NHC content in bulk [EMIM+][OAc–] to be on the order of parts-per-million (ppm) levels, with order-of-magnitude enhancement of NHC concentration at the liquid/vapor interface. The interfacial enhancement of NHC content is due to both poorer solvation of the ionic reactants and solvophobic stabilization of the neutral NHC molecule at the liquid/vapor interface.


Procedure for generating [AcOH][OAc − ] dimer training data
As mentioned in the main text, we run AIMD simulations with CP2k S1 to generate training data for proton transfer between acetic acid and acetate. The QUICKSTEP method was used, with the PBE-D3(BJ) exchange-correlation functional and a mixed Gaussian/planewave basis set. S2-S4 The aug-TZVP basis set, a plane-wave cutoff of 280 Rydberg to expand the auxiliary electron density and Goedecker-Teter-Hutter pseudopotentials for the core electrons were used. S3 SCF convergence criteria was set to 10 −7 au. NVT simulations were run at temperatures of 300K, 400K, 500K and 600K using a Langevin thermostat with a 0.005 fs −1 friction coefficient and a 1 fs timestep. As the proton is shared closely between the two acetate molecules during the gas-phase AIMD simulations, three additional simulations were run at 300K with harmonic umbrella potentials. This allows for sampling dimer geometries at greater intermolecular separations. The umbrella potentials were placed between an oxygen atom on one molecule and the reactive proton. Each umbrella potential ( 1 2 kx 2 ) used a 200.0 kJ/mol/Å 2 force constant and were centered at values of 3, 3.25 and 3.5Å.
The AIMD simulations yielded approximately 290,000 acetic acid/acetate dimers. We computed SAPT0/aug-cc-pVTZ energies for individual configurations with Psi4. S5 As this is a symmetric reaction, two SAPT0 calculations were performed for each geometry: one with the acetic acid classified as the molecule with the shorter minimum O-H bond distance and one with the acetic acid classified as the molecule with the longer minimum O-H bond distance. This yields 580,000 SAPT0 calculations. These were used to train the E Inter,N N Solute neural network. We use AP-Net as in our prior work for the neural network architecture. S6,S7 Details of the AP-Net neural network hyperparameters and training details are described in the next section. As the E Inter,N N Solute term is a residual correction to E N onbonded,F F Solute,Solvent , we plot E Inter,N N Solute + E N onbonded,F F Solute,Solvent vs. SAPT0δ HF energies for a subset of the test set in Figure  S1. The mean absolute error (MAE) for the total test set is 0.55 kJ mol −1 . As mentioned S4 previously, E Inter,N N Solute is multiplied by a damping function. We use the Fermi-Dirac function for this term; the parameters are shown in Table S1. S5 2 E Inter,N N Solute

AP-Net hyperparameters
The AP-Net architecture is used for modeling the E Inter,N N Solute term. As in our prior work, S6 the interaction energy predicted by AP-Net is modified slightly from the published definition by Glick et al. S7 The pairwise energy and sum of pairwise energies are multiplied by separate damping functions, as shown in Equations S1-S3: S2) f (r bond ) = 1 (e β(r bond −µro) + 1) (S3) The AP-Net sum runs over atoms on monomers A and B. In Equation S1, the f c (r ij ) is a standard cutoff function used in many neural network frameworks. S8 The r c cutoff is set to 4.0Å, which ensures that there is no contribution to the total interaction energy by atom pairs separated by more than this amount. This restricts E Inter,N N Solute to serve as a short-range correction. The other cutoff function, f (r bond ), is a Fermi-Dirac function. The function takes as input the length of the dissociating bond. If the bond stretches too far, then the contribution to the total energy from the diabat is small. The Fermi-Dirac function zeros out the energy of the neural network at these distances. The parameters for the cutoff function are shown in Table S1. The r o term is listed as the equilibrium bond length for the bond in question, and µ is an arbitrary factor that scales it by some amount. The procedure for selecting µ involves determining where the difference between the diabat and reference energy becomes zero, and is explained further in the Supporting Information of Stoppelman and McDaniel S6 network from SchNet was used as the activation function. S9 The loss function is the mean squared error (MSE) loss function shown in Equation S4. The loss is minimized between the (SAPT0δ HF ) -SAPT-FF energies that make up the training data and the neural network prediction.
For training the neural network, we use a train-validation-test split of 80:15:5. The neural network is trained using the Adam optimizer, with an initial learning rate of 5.0 x 10 −4 and a decay rate of 0.8. S10 Training finished once the learning rate decayed to a value of 1.0 x 10 −6 . We use batch sizes of 100 for both training and validation. The mean absolute error (MAE) for E Inter,N N Solute on the test set is shown in Figure S1.

H ij AP-Net hyperparameters
The modified version of AP-Net still models pairwise energies between atoms, but the restriction of the pairs being on opposite monomers is removed.
The  S6) ρ is set to 0.1 here.

Configuration
For condensed phase simulations, we found it necessary to modify the E Intra,N N Solute acetic acid term from its original parameterization. S6 In the gas-phase, the syn configuration of acetic acid is the stable form; however, in water, the syn configuration has been found to be approximately isoenergetic with the anti configuration. S11 We find that this is true for the bulk ionic liquid as well. As our neural networks were built with structures mainly taken from gas-phase AIMD simulations, the anti configuration was not properly represented in each neural network's training set. In order to add these structures to the various training sets, we ran gas-phase AIMD simulations of the NHC/acetic acid dimer and the acetic acidacetate dimer with the acetic acid dihedral restrained to various values from 0 • to 180 • .
We took approximately 50,000 acetic acid configurations from these simulations, computed PBE-D3(BJ)/aug-cc-pVTZ energies/forces and added them to the E Intra,N N Solute acetic acid training set; training details for this can be found in the prior work. S6 SAPT0/aug-cc-pVTZ calculations were then run for both NHC/acetic acid and acetic acid/acetate dimers (again approximately 50,000 structures for both) and then both E Inter,N N Solute neural networks were retrained. Finally, H 12 and H 23 were both retrained with DFT energies and forces on the collected dimers. We find that this is sufficient for modeling the varying acetic acid dihedral angle configurations in the reaction.

S10
5 Acetic Acid -Acetate Umbrella Sampling Procedure We employ umbrella sampling to produce a 2D free energy surface for proton transfer between acetic acid and acetate. The collective variable we used is shown in Equation S7: This is the minimum distance between the oxygens on one molecule with the reactive proton. This CV is computed for both molecules in order to construct the free energy surface.
Umbrella sampling simulations were run for both AIMD and PB/NN using the Hamiltonian described in Section 2.1 of the main text. The AIMD simulations were run using the CP2k/ASE interface, S1,S12 and the PB/NN simulations were also run with ASE. Plumed was used to apply the harmonic umbrella potentials. S13 A 200.0 kJ/mol/Å 2 force constant was used for each umbrella ( 1 2 kx 2 ). The umbrella potentials were centered at values shown in Table S3. Each window was run for approximately 20 ps. The free energy surfaces are shown in Figure 2 in the Main Text.
In the main text, we discussed the issue of "reinitializing" the PB/NN Hamiltonian when the reactive complex identity changes during the course of the simulation. This occurs when a solvent acetate ion "OAc − s " displaces a closer contact acetate ion "OAc − 2 " in the reacting complex. Here, we provide more detail for how the PB/NN Hamiltonian is reinitialized in this case. ); this would need to be shortened to 2 A or less in order to avoid having the "switching" procedure affect the total ground-state energy. We find that this would hinder the ability of the E Inter,N N Solute models to function as a correction to the E N onbonded,F F Solute,Solvent , so we do not adjust the cutoff. The fluctuations in the total energy caused by the switching procedure are on the order of 10 kJ mol −1 and can be observed in Figures S4d and S5d.
The last set of terms that would need to be adjusted is the coupling terms. In this scenario, the H 13 and H 23 coupling terms would need to be switched from OAc − 2 to OAc − s . Again, since the separation between OAc − 2 and the reacting proton is large, the value of both coupling elements will be close to zero and thus switching them between molecules will have no effect on the total energy. S14

PB/NN Energy Conservation
Here we show the energy conservation of the PB/NN Hamiltonian, comparing it to a standard simulation. Figure S3 shows the extent of energy conservation for a nonreactive NVE simulation of 40 [EMIM + ][OAc − ] pairs. Standard force field terms are used for the intramolecular components, and the SAPT-FF force field is used for the intermolecular terms. The simulation is run using ASE calling OpenMM as a calculator. Due to the Drude SCF procedure, there is an energy drift of approximately 16 kJ mol −1 ps −1 . In Figures S4 and S5, we show the PB/NN energy conservation for different points along the reaction profile shown in Figure 2 of the main text. Figure S4a and Figure S4b show the energy conservation for windows with umbrella potential values of ∆(CH − OH min ) = −1.0Å and ∆(OH min,1 − OH min,2 ) = 0.0 A. Figure S4c and Figure S4d show the energy conservation of windows with umbrella potential values of ∆(CH − OH min ) = −1.0Å and ∆(OH min,1 − OH min,2 ) = 2.5Å. Note that Figure S4a-b essentially show the same extent of energy conservation as Figure S3. This is because the proton is localized on EMIM + , and thus the ground-state energy is equal to H 11 . Also, the umbrella potential fixes the acetates to be close to equidistant to the reactive proton, meaning that the switching procedure does not need to be employed for this window. Figure S4c-d show various jumps in the total energy due to the switching procedure for E Inter,N N EM IM + /OAc − as acetates from the solvent may come closer than the acetate in the reacting complex. Figure S5a-b show the energy conservation for a simulation with umbrella potential values of ∆(CH − OH min ) = 3.5Å and ∆(OH min,1 − OH min,2 ) = 0.0Å.
Note for these windows the ground-state energy will be a mix of H 22 and H 33 , as the proton is shared equally between an acetic acid molecule and acetate. This leads to slightly worse energy conservation than seen in Figure S3 and Figure S4a-b. Figure S5c and ∆(OH min,1 − OH min,2 ) = 2.5Å. The switching procedure is applied again here, as it is possible for solvent acetates to approach closer than the acetate considered part of the reacting complex. The jumps in energy are smaller in this case, as now the E Inter,N N AcOH/OAc − S15 model is switched between acetate molecules. These are essentially limiting cases of the energy conservation performance; Figure S4a-b represent the best possible performance for the Hamiltonian for windows where the acetates are placed equidistant to the proton. Figure   S4c-d places the secondary acetate as far away from the proton as is possible for our reaction profile, so this would lead to the largest number of events that requires switching acetates in and out of the Hamiltonian. Every other window would exhibit energy conservation intermediate to these two simulations. All these fluctuations are smaller than those caused by the NVT thermostat, but we do acknowledge that the procedure may cause some (relatively minor) artifacts. For c) and d), this umbrella potential is set to 2.5Å. Figure S5: Energy conservation for umbrella sampling potentials of ∆(CH − OH min ) = 3.5 A. The secondary umbrella potential, ∆(OH min,1 − OH min,2 ), is set to 0.0Å for a) and b). For c) and d), this umbrella potential is set to 2.5Å.

Diabatization Limitations
As we noted in our previous work, the current scheme for the diabats means that the intermolecular interactions for the reacting complex change from a DFT description to SAPT0 δ HF description as the reaction occurs. The region where the switch between the two levels of theory happens is determined by the Fermi-Dirac functions applied to the H ij terms.
There is some ambiguity in determining the center for the Fermi-Dirac function, although we aim for the point along the reaction coordinate for which the difference between DFT and SAPT0δ HF interaction energies is close to zero. Unfortunately, there is no guarantee that DFT and SAPT0δ HF will agree for every set of configurations or that they will agree after the difference becomes zero. We show potential energy scans in Figure  The geometry optimizations were performed using the geomeTRIC software package interfaced with Psi4. S5,S14 In Figure S6 a EMIM + /acetate dimer is modeled; in Figure S7 a NHC/acetic acid dimer is modeled; and in Figure S8 an acetic acid/acetate dimer is modeled. Asymptotically, the energies eventually agree with one another, but there is a significant discrepancy at intermediate distances. One could increase the range where H ij is applied until DFT and SAPT0δ HF interaction energies agree asymptotically, but this would delay switching to the pure diabatic state. Physically reasonable diabat definitions are important for modeling solvent effects on chemical reactions in these methods. The differing levels of theory do cause a sudden increase in the potential energy at the transition to the pure diabat, which can be seen in Figures 3, 4 and 9 in the Results section and also Figures S6-S8; although the transition is smooth due to the damping functions, this represents an unphysical artifact on the simulations. As we are focused on environmental effects on reactions 1 and 2 in this work, this artifact fortunately does not affect this analysis, as it occurs regardless S19 of environment; however, it does impact the quantitative results. We treat our results as semi-quantitative, and leave determining a more rigorous/effective diabatization procedure for further work. We do slightly extend the Fermi-Dirac damping function for E Inter,N N Solute between EMIM + /acetate and the Fermi-Dirac damping function for H 12 that takes the EMIM + C-H bond distance as input (see Equation S5 for the form of the coupling terms).
This extends the point along the reaction for which the Hamiltonian transitions to the pure H 22 state, which slightly decreases the change in potential energy at the transition, as shown in Figure S8. The new vs. old damping parameters are shown in Table S4. The energies are plotted with respect to the distance from the closest acetate oxygen to the reactive proton. Figure S7: Comparison of the total DFT energy, monomer DFT energy + SAPT0 energy, monomer DFT energy + SAPT0δ HF energy and PB/NN energy for NHC/AcOH. Note the difference between the relaxed, isolated gas-phase monomer energies of NHC/Acetic acid and EMIM + /acetate is added to the PBE-D3(BJ) monomer + SAPT0 and PBE-D3(BJ) monomer + SAPT0δ HF terms. The energies are plotted with respect to NHC carbonproton distance. Figure S8: Comparison of the total DFT energy, monomer DFT energy + SAPT0 energy, monomer DFT energy + SAPT0δ HF energy and PB/NN energy for AcOH/OAc − . The energies are plotted with respect to the distance from the closest acetate oxygen to the reactive proton.

Gibbs Dividing Surface
We plot the number density profile along the z-dimension for a 100 ns liquid/vapor interface (nonreactive) simulation, using the OpenMM settings described in the main text. We count the number of ions (separately for cations and anions) within discretized bins along the zaxis; we use the center of mass for each ion belonging to each type for the binning procedure.
We divide the simulation cell in two and average the results from the center of the cell. Note that the liquid is placed at the center of the unitcell and there is a vacuum gap on both sides.
The dimensions of the cell are 22.08Å x 22.08Åx 66.24Å. The Gibbs dividing surface occurs where the density falls to approximately half that of the bulk value, which is approximately 11Å from the center of the simulation cell in this case. Figure S9: Number density profile for both EMIM + and acetate at the air-liquid interface.

S22
10 Umbrella potentials used for gas phase dimer Umbrella sampling for the gas phase dimer was performed using r CH and r OH min (the minimum oxygen-proton distance) as the collective variables; the center umbrella sampling windows are listed in Table S5 and a force constant of 200.0 kJ/mol/Å 2 was used for all umbrella potentials ( 1 2 kx 2 ). These simulations were run with a 2x2 Hamiltonian (upper 2x2 block of Equation 3 in the Main Text), as there is no third reacting species for the dimer.  Figure S10 shows the liquid phase PMF with labels for different configurations along the profile. Reactive complex trimer geometries are shown in Figures S11, S12 and S13 for the liquid, gas phase and ionic liquid/vapor interface. Note that the liquid and gas phase configurations are noticeably different from one another past point a) of the profile, while the liquid and liquid/vapor interface configurations are fairly similar. For example, in Figure   S13c the nonreactive spectator acetate for Reaction 1 in the reactive complex is observed to have a similar proximity to the EMIM + ethyl group as we see in the liquid phase simulations Figure S11c. This acetate was positioned over the imidazolium ring in the gas phase, as can be seen in Figure S12c. S24 Figure S11: Simulation frames showing proton transfer from EMIM + to the AcOH/OAc − dimer in the bulk liquid. The location of each panel a-f on the reaction profile can be seen in Figure S10.
S25 Figure S12: Simulation frames showing proton transfer from EMIM + to the AcOH/OAc − dimer in the gas phase. The location of each panel a-f on the reaction profile is the same as in the liquid phase and can be seen in Figure S10. Figure S13: Simulation frames showing proton transfer from EMIM + to the AcOH/OAc − dimer at the liquid/vapor interface. The location of each panel a-f on the reaction profile is the same as in the liquid phase and can be seen in Figure S10.

S27
12 NHC ring center of mass / AcOH -OAc − Dimer center of mass distance The histograms in the left column of Figure S14 are of the distance between the NHC center of mass and the AcOH/OAc − dimer center of mass from the liquid and gas phase. The histograms in the right column are of the angle between the C2 carbon, the NHC center of mass and the AcOH/OAc − dimer center of mass. The collective variable values that each histogram was constructed from (with the collective variables defined in Equations 5a and 5b of the Main Text) are listed in each figure title. Figure S15 shows the same plots but for the liquid phase and liquid/vapor interface. The histograms are noticeably different for the gas and liquid phase, while they are similar for the liquid phase and the liquid/vapor interface. Note that the angle distributions are broader in both the gas and liquid/vapor interface than in the liquid phase. This due to the absence of solvating ions in the gas phase and lower ion density at the interface. Figure S14: Each row contains two histograms per window. The first column is a histogram of the NHC center of mass and AcOH/OAc − dimer center of mass distance. The second histogram in the row is a histogram of the angle between the C2 carbon, the NHC center of mass and the AcOH/OAc − center of mass. Results are shown for the liquid and gas phase. Figure S15: Each row contains two histograms per window. The first column is a histogram of the NHC center of mass and AcOH/OAc − dimer center of mass distance. The second histogram in each row is a histogram of the angle between the C2 carbon, the NHC center of mass and the AcOH/OAc − center of mass. Results are shown for the liquid and liquid/vapor interface.

Spatial Distribution Functions
We plot spatial distribution functions of the EMIM + ring surrounded by oxygen atoms for the collective variable values of (CV 1 = −1.0Å, CV 2 = ±2.0Å); the SDFs were constructed using Travis. S15 The different configurations configurations for these two windows leads to the asymmetry in the liquid PMF.
We also show SDFs from the liquid and liquid/vapor interface in Figure 11 of the Main Text, with both corresponding to the R 1 configuration. Except for the nonreactive ring proton near the ethyl group, the oxygen density is similar.

S33
14 EMIM + H3 and H4 ring protons/oxygen RDF Figure S17 shows ρg(r) between the (reacting complex) EMIM + nonreactive H3 and H4 protons and the acetate oxygens for 3 different umbrella sampling windows in the bulk liquid. The window constraints are listed in the Figure S17 caption. 15 Ring proton -oxygen RDFs In Figure S18, we show the combined H3/H4 ring proton/oxygen atom ρg(r) distributions for both the "Reactant" EMIM + (R 1 ) and "Product" NHC (P 2 ) states mentioned in the Main Text. The hydrogen bond between these molecules is much stronger for the Reactant.
In Figure S19, we show the Reactant EMIM + and Product NHC ρg(r) decomposed between H3 and H4. Figure S19b shows that there is small acetate presence at the H3 proton and negligible presence at the H4 proton for the NHC. These simulations are taken from the liquid phase.  Figure S18: ρg(r) for the reactive EMIM + H3 and H4 proton/acetate oxygen atoms and ρg(r) for the NHC H3 and H4 protons/oxygens. The "Reactant" simulation is from a window with (CV 1 = −1.0Å, CV 2 = 0.0Å) and the "Product" simulation is from a window with (CV 1 = 3.5Å, CV 2 = 0.0Å).  Figure S19: a) RDF of the reactive complex EMIM + H3 and H4 protons with solvent acetate oxygen atoms. This is a replot of Figure S17c. b) RDF of the reactive complex NHC product H3 and H4 protons with solvent acetate oxygen atoms.

Reacting Complex Atom Labeling
We follow the atom labeling convention in the reactive complex for some sections of the paper. The reactive complex is considered to be composed of EMIM + and two acetates, and allows for solvent acetates to switch in and out of the reactive complex. Figure S20: Reactive complex and atom naming convention in this work.

S37
In order to assess the relative stability of the NHC at the liquid/vapor interface in comparison to the bulk liquid, we performed nonreactive umbrella sampling simulations of the NHC in the liquid/vapor interface system. The umbrella sampling simulations are of a system consisting of the NHC, acetic acid and 39 [EMIM + ][OAc − ] ion pairs, and are nonreactive; this prevents simulation of the AcOH/OAc − dimer, but allows for focusing on the environmental preference of the NHC. The z-coordinate of the NHC center of mass was biased from the center of mass of the simulation box (set at the origin, 0.0Å) to 20Å in 0.5Å increments. A 20 kJ/mol force constant was used for the harmonic umbrella potentials. The simulation settings are the same as those mentioned in the main text. The 1D PMF is shown in Figure S21. The minimum free energy region is at ∼ 11Å, which is approximately the location of the Gibbs dividing surface. There is a 15-20 kJ/mol barrier for the NHC to move further into the liquid or into the vacuum.
From our results, it appears that there is a positive solvation energy for the NHC in the IL. However, it has been found experimentally that aromatic compounds such as benzene are solvated to a limited extent in [EMIM + ][OAc − ], which implies a small solvation energy for the NHC here. S16 Determining the solvation energy of NHC within statistical accuracy would require more extensive simulation, which is beyond the scope of this work. Figure S21: Potential of mean force for the NHC center of mass along the z-dimension of the simulation box.

Estimate of NHC Concentration in the Bulk Liquid
Here we give an order of magnitude estimate for the NHC concentration based off of our free energy results for the liquid and liquid-interface systems. The equilibrium coefficient is given by the ratio of the forward and reverse rate constants for the reaction: Considering the reaction: the forward and reverse rates are given through second-order kinetics as: with the forward and reverse rates equal to one another at equilibrium. Additionally, the S39 rate constants can be related to the free energy through the well-known expression: The ratio of the forward and reverse rate constants allows for the relation of the concentrations of the species to the free energy: where we have assumed that the prefactors in Equation S15 and S16 are approximately equal to one another. The concentration of EMIM + and OAc − in the bulk IL is equal to 6.1 M.
Entering this in to Equation S17 and equating the concentration of NHC with that of acetic acid gives: [N HC] ≈ (e −70/2.5 6.1 2 ) 1/2 = 5.1µM (S18) in the liquid phase and [N HC] ≈ (e −60/2.5 6.1 2 ) 1/2 = 37.9µM (S19) for the liquid/vapor interface. Note that the free energies above are approximate and read off of Figure 9 in the main text. These concentrations should be treated as order of magnitude estimates rather than exact figures. is the slowest part of the simulation and will scale with respect to the size of the system. This is due to various inefficiencies in the current code implementation associated with sending the full system positions and forces back and forth from the CPU and GPU.

PB/NN Timings
A 1 ps simulation of the bulk liquid system described in the Main Text using the PB/NN Hamiltonian takes 833.7 seconds to complete. For comparison, we also ran a 1 ps simulation of a nonreactive polarizable force field simulation of 40 [EMIM + ][OAc − ] pairs (essentially E Bonded,F F Solvent + E N onbonded,F F Solute,Solvent ). OpenMM was used for this simulation. S17 The Drude SCF Integrator was used for this system with a timestep of 0.5 fs; an anharmonic restraining potential was used for the Drude oscillators. S18 The OpenCL OpenMM kernels with mixed precision were used for this simulation. Other simulation settings are the same as those used for the reactive simulations (mentioned in the Methods section of the Main Text). This simulation ran for 141.9 seconds, which is 5.88 times faster than the reactive PB/NN simulation. A code requiring a lower extent of CPU/GPU communication would reduce the differential between the two simulations, and we will pursue this in future work. For comparison, we also ran simulations with the nonreactive polarizable force field using the Drude Langevin integrator, S41 also run using the OpenCL kernels in OpenMM. S18 This simulation was approximately 180 times faster than the PB/NN simulation. Note that the PB/NN procedure uses the Drude SCF method to calculate the polarization energy, which is less computationally efficient than the extended Lagrangian method used in the Drude Langevin integrator.